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f"**) Abstract 

Quantum computers could potentially simulate the dynamics of systems such as 
polyatomic molecules on a much larger scale than classical computers. We investigate 
i a general quantum computational algorithm that simulates the time evolution of an 

arbitrary non-relativistic, Coulombic many-body system in three dimensions, consider- 
ing only spatial degrees of freedom. We use a simple discretized model of Schrodinger 
evolution and discuss detailed constructions of the operators necessary to realize the 
scheme of Wiesner and Zalka. The algorithm is simulated numerically for small test 
i— i cases, and its outputs are found to be in good agreement with analytical solutions. 

+1 1 Introduction 

g 

Quantum computers could provide an exponential speedup over classical methods of sim- 
ulating quantum mechanical systems pQ. Wiesner [2] and Zalka |3j have sketched a basic 
theory of such simulations that includes discretization techniques and decomposition of time 
evolution operators. Other approaches in the literature include the quantum lattice gas au- 
tomaton (QLGA) |11 | I12 | [13] and various methods of decomposing fermionic Hamiltonians 

\Q into spin operators via the Jordan- Wigner transformation [HJ [141, \TE[ [T6] . 

However, work based on the first approach has generally focused on bounds for relevant 
time or space complexities rather than the explicit form of the requisite operators [U [T7] , 
as well as on the problem of calculating energy spectra [16] . Conversely, in this note, we 
expand on Wiesner and Zalka's approach to compute position-dependent wavefunctions. 
We propose detailed constructions of the operators (quantum gates) for such simulations 
and numerical tests of the algorithm described. 

ITpi Given the initial state of some number of particles confined to a finite volume in three- 

dimensional space, we aim to compute the wavefunction of the system after an arbitrary 
period of time. Our simplified model uses only the non-relativistic Coulomb Hamiltonian, 
whose terms include kinetic energies and electric potentials, and neglects spin-orbit effects. 



2 Method 
2.1 Hamiltonian 

The Hamiltonian H is a Hermitian operator acting on states in a Hilbert space. For 
quantum computation, the Hilbert space in question is finite-dimensional. We will consider 
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primarily the Coulomb Hamiltonian H = T + U , which consists of terms of the form 

T — — V 2 and U — ^ 

' 2—* 2rrii 1 ~ 47re ^ ^ ||r; - r,-|| ' 

i i j>i J 

For a molecular system, 

H = f e + f n + U ee + U nn + U en (1) 

where the U terms represent various interactions between electrons and nuclei. One of 
our goals is to make such continuous operators precise in the context of quantum comput- 
ing: what does the single-particle momentum operator —ihV mean as a finite-dimensional 
unitary matrix and a quantum gate? 

We use the following simulation parameters to determine the dimension of the Hilbert 
space and the evolution of ip: T, the total time over which to simulate; Nt, the number of 
small time increments; L, the bound on the dimensions of the system; and n, a parameter 
chosen to reflect our desired precision. 



2.2 Discretization 

Wiesner [2] and Zalka [3] have outlined algorithms for simulating quantum many-body 
systems with the common element of diagonalization using the quantum Fourier transform 
(QFT). Expanding on their framework, we explicitly construct the operators needed to 
carry out these simulations. 

We first review the common encoding scheme for Cartesian coordinates in three spatial 
dimensions. Let all particles be confined to a finite volume (0 < x, y, z < L). We divide the 
intervals along each coordinate axis into 2 n subintervals of length 5 = L/2 n . The volume 
is hence partitioned into many subvolumes (position basis states) where each amplitude 
of the state vector of a particle in the system specifies the probability that it will be 
found in a subvolume. The system wavefunction is a tensor product of discretized single- 
particle wavefunctions, or a sum thereof, encoded in a register of qubits. We approximate 
a continuous single-particle wavefunction i[)(r,t) by the state vector 

2 n — 1 2 n -l 2"-l 

|V(r,i)) = E E E H^k,t)\ijk), (2) 

i=0 j=0 k=0 

normalized by a factor of j k l^( r j?fc> ^)| 2 ) -1 ^ 2 ; where 

njk = (r x ,r y ,r z ) = 5 (i + 1/2, j + 1/2, fc + 1/2) . (3) 

The one-dimensional position basis states and \k) range over discrete positions on 

the x, y, and z axes, respectively, while \ijk) denotes the tensored state \i)\j)\k). We assign 
the amplitude associated with each subvolume to be the amplitude of tp at its center. Each 
particle requires a register of n qubits for each of the x, y, and z intervals. The algorithm 
thus requires 3nN total qubits, N being the number of particles. 

The particles are assumed spinless, but that leaves the question of distinguishability. 
Our algorithm is fully compatible with simulating indistinguishable fermions if one adds 
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ancillary qubit registers to allow for the initialization of an antisymmetric spatial wavefunc- 
tion via the algorithm of Abrams and Lloyd [8]. In our notation, this requires 0((3nN) 2 ) 
operations given N single-particle states. This procedure also suffices (and requires fewer 
operations) for symmetrization — for example, when dealing with bosonic nuclei. Our main 
focus here is the form of the operators, so we will not pursue this point further. 
From the Schrodinger equation Hip = ihdtip, we have the time evolution 

|^r,t + e)> = e-^ e |^Mo)>. (4) 

For very small time scales e, we can apply the Trotter decomposition |17j to construct gates 
for each term of the Hamiltonian separately and apply them sequentially as unitary time 
evolution operators: 

e -s^ = e -s^ e -^ + 0( e 2). (5 ) 

To carry out the algorithm, we choose a time interval T over which to simulate and partition 
it into Nt small intervals e = T/Nt- At each of the Nt time steps, we apply the operator 
([5]) to the state vector of our system in the position basis. 

Naively, one can repeat the entire algorithm many times and make subsequent mea- 
surements of every qubit to iteratively determine a probability distribution for the particle 
configurations, thus allowing one to predict the overall spatial density at the end of the 
desired time interval. The procedure is identical to that required for the algorithm of 
Abrams and Lloyd [8] for simulating fermions in the first-quantized representation, namely 
constructing a histogram of configurations obtained from the repeated measurements. 



2.3 Construction of Kinetic Energy Operators 

It is commonly noted [11 [3l SI [T7] that implementing the kinetic energy operator by 
Fourier transforming to the momentum basis and back reduces the time evolution operators 
to diagonal operators and the QFT. However, the form of the kinetic energy operator in the 
model of Wiesner and Zalka has not, to the author's knowledge, been previously addressed. 

To this end, we determine the form of the momentum operator matrix P by analogy 
with the continuous case. First, consider the —iftV operator for a single particle in one 
dimension: i.e., along an x-interval of length L divided into 2 n subintervals of length 5. 
To obtain approximations for the derivative of the discrete-valued wavefunction ip at each 
point k (0 < k < 2 n ) corresponding to a position basis state, define 

n+ _ j,(k + l)-^(k) n ave_ ^(fc + l)-V>(fc-l) 

° k ~ 6 ' k ~ 6 ' k ~ 25 

so that to the best approximation, 

Ad^(O) 
D x ^(2 n - 1) 
D x il>(k) 

The action of the "discrete" derivative operator is thus described by 

2™— 1 „ 2 n — 2 

£ ip(k)\k) ^ D+10) + £ Dr\k)+D^\2 n - 1), 

fc=0 k=l 



= D+ + 0(5), 

= D2n_i + 0{5), and 

= Dt ve + 0(5 2 ) (0<A:<2 n -l). 
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or in matrix form as 
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(6) 



In addition, we require that the matrix — |(P 2 /2M)e be skew-Hermitian so that its expo- 
nential is unitary. Hence P 2 must be symmetric, and we can make it so by sacrificing the 
approximations for D x at the endpoints. The corresponding momentum operator is 



P = —ihD x ~ — 



ih 
26 



( o 



1 




••• 0\ 



-i oj 



(7) 



The form of the kinetic energy operator follows immediately, and extending the kinetic 
energy operator to three dimensions is as simple as applying it to the qubits encoding the 
x-, y-, and z-coordinates of each particle separately. Since this construction takes place 
in the position basis, its implementation does not require a QFT. We next discuss two 
methods of implementing the kinetic energy term of the time evolution operator e~^ Te as 
a set of quantum gates. 

As a first method, note that 




■|0)<0|+ A-|2"-l)(2 n -l| 



i=i 



where P; 



+ 1| - 2|i)(i| + \i + - 1|. Let £ = The Trotter formula gives 



\2M 



-C|0>(0| 



TJ e ^ j e -€|2»-i><2»-i| + 0( £ 2) 



(8) 



where e -«l ><°l and e -€|2"-i><2»-i| 

are both single-qubit phase rotation operators. Further- 
more, since each of the Pi is a very sparse matrix, the e^ Pi are block diagonal matrices 
representing two-qubit operations: e^ Pi = diag(/j_i, Mp, Ii n -2-i) where 

/ i 

M P = exp -2£ 

V e o o 

and I m denotes the m x m identity matrix. Hence we have decomposed the kinetic energy 
operator entirely into controlled one- and two-qubit gates. Since 2 n such gates act on each 
n-qubit register representing position basis states along one coordinate axis, 3N2 n one- 
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and two-qubit gates are needed to simulate the kinetic time evolution operator of the entire 
system. This gate complexity is linear in the number of subintervals along each axis, or the 
"absolute" precision N = 2". 

A second, and more robust, method is to realize that in the continuum limit, P will 
be exactly diagonalized by the QFT, as follows: letting D be the dimension of the Hilbert 
space, F be the QFT of dimension D, and a = —ih/28, we can write 

D-2 



P = a ^2(\x)(x+ 1| - \x + l)(x\). 



1=0 



Applying F = -i= Yljk=o e- l27T -' k ^ N \j}(k\ and using (j\k) = 5jk for simplification shows that 



fN ^0 

D- : D- ! 

D 



Fpp] = ji E T,^ klDs ^M-^ k/D S2{.jM\3){M (9) 



j=0 k=0 

where 

D-2 D-l 

S ± {j, k) = J2 e i27rp( - j ~ k ^ D and S 2 {j, k) = e i2 ^~ k ^ D . 

p=0 p=l 

Summing the geometric series gives Si(j, k) = —e i27r v—k){D—i)/D anc j g 2 (J f k) = —1 if j ^ k. 
Consequently, as D — > oo, S\(j, k) — > —1 and the non-diagonal entries approach 

a^k/D_ e -i2*k/D )= *L isin ^\^ Q (1Q) 

If j = k, Si(j, k) = Sz(j, k) = D — 1, so that as D — > oo, the k th diagonal entry becomes 

D ^a(e-^ k l D - e l2 ^ D ) -> -2iasin — = -^sin— . (11) 
D ' D 5 D y ' 

Thus FPF^ is a diagonal matrix with real eigenvalues in the continuum limit, as expected 
physically. Now that we have explicitly determined these eigenvalues, we can construct the 
kinetic energy operator from QFTs and a diagonal P matrix with entries of the form (11) in 
the momentum basis. Once again, an unoptimized method of implementing the associated 
time evolution operators requires only 3iV(2 n + 0(n 2 )) gates — taking into account both 
the QFT and the diagonal component — which is linear in both the number of particles 
and the number of subintervals along each coordinate axis. 



2.4 Construction of Potential Energy Operators 

The construction of the potential energy operators is generally straightforward in the po- 
sition basis. For example, the Coulomb potential operator for a molecular system has the 
following physical parameters: N e ,ri (resp. iV n ,Rj) for the number and positions of the 
electrons (resp. nuclei) and Zi,M$ for the atomic numbers and masses of the nuclei. Take 
r p - r i P jpk P (resp. R P = R/pjpft-J as shorthand, where i p j p k p (resp. I p J p K p ) is a bit string 
corresponding to the binary representation of some index describing a position basis state, 
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or subvolume, occupied by the p th electron (resp. nucleus). Recall that the position vector 
rjjfc is defined in equation ([3]). Also, let A[m] denote the m th diagonal entry of matrix A. 

U ee can be seen as an operator on the state of the entire system of electrons, or as a 
23nN e x 2^ nN e diagonal matrix whose entries give the total potential energy due to any of 
the 2 3nNe allowed configurations of electrons. If 

\ m ) = \isj 8 k 8 ) 

5=1 

represents one of the many possible electron configurations, then 

p=l q=p+l r H 

Analogously, U nn acts on the state vector of the entire system of nuclei and U en on all 
particles. For those diagonal entries that represent configurations of the system in which 
multiple particles are found in the same "box" of volume (i.e., when \i p j p k p ) = \i q j q k q ) 
and hence r p = r q for some p, q), we can set the potential to some approximate large, but 
finite, value. For example, in the case of two electrons, we can replace the undefined term 
e 2 /47reo||r p — r g || with the approximation e 2 /47reod. 

2.5 Complexity of Potential Energy Operators 

It is unknown whether nonlocal potential energy operators such as these can be implemented 
in a polynomial number of resources; Lloyd demonstrated efficiently-scaling implementa- 
tions for local interactions [7|. Benenti and Strini [1] have shown that the time evolution 
operators corresponding to a harmonic oscillator potential do scale polynomially with the 
number of qubits, and their argument is easily extended to power-law potentials of the 
form V(x) = ax r where r is a positive integer. The difficulty with the Coulomb potential 
lies in r being negative, which prevents easy factorization of the operators into products of 
phase-shift gates acting on a fixed number of qubits, even after introducing ancillary qubits. 
However, it may be possible to implement these nonlocal potential operators in polynomial 
time if we choose to sacrifice some additional precision in computing the values of our po- 
tential energies. Polynomial scaling is important because it is necessary to first calculate 



the potential energies ( 12 ) on a classical computer before implementing the corresponding 
phase rotations in quantum gates. 

First, we can utilize redundancy to implement the potential energy operators more 
efficiently. As an example, consider the three-qubit phase shift operator 

diag(e iei , e id2 , e iH , e i9i , e Wa , e id4 , e if>1 , e id2 ) 

with "redundant" entries. The naive method of constructing this diagonal operator is to 
use the first two qubits as controls for phase shift operators on the third qubit, resulting in 
four gates acting on the third qubit. A more efficient scheme is shown in Figure [TJ using 
the single-qubit gates 

a ( eWl ^ a a ( ^ 
012 = { o & and 6u = [ e *« 
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Because there are many configurations of a molecular system with the same potential energy 
(up to a precision), the potential energy operators also have many redundant entries. In fact, 
let U be one of U ee , U nn , and U en , and let the corresponding N be N e , N n , or N e +N n . Then 
U[x] = U[2 3nN — 1 — x] for < x < 2 3nN ~ 1 — 1, meaning that the potential energy operators 
are symmetric about the antidiagonal. This is easily seen by noting that the indices of the 
states \x) and \2 3nN - 1 — x) are bitwise complements, meaning that by symmetry, they 
represent reflected versions of the same molecular configuration. This observation reduces 
the number of gates needed by a factor of two, and similar optimizations could lead to 
much more efficient implementations of the potential energy operators. 

\ x i) ~T T~~ 



\x 2 ) -e- 



\ X V 



e — t — CD cb 



#34 



H2 



Figure 1: Efficient implementation of a diagonal operator. 

Second, because we would like to simulate arbitrarily large systems, the following vari- 
ables are expensive: N e , N n (in addition to the atomic numbers Zi), and the absolute 
precision ./V where n = log 2 N. We have already seen that the kinetic energy operators 
scale linearly with N. The number of distinct values of the potential energy of the system 
can scale polynomially if we restrict our calculations to some appropriate precision AU. 
Then, despite the exponential number of particle configurations, their potential energies 
only assume a polynomial number of values, perhaps allowing us to employ redundancy to 
implement the potential energy operators efficiently. 

To see how this might work, let e' = e 2 /47reo- The extreme values of the potential 
energy of the entire system might go like 



^^l^+EE^I andC/ min = -U>. 

j>i 



s 



We would like to determine some increment AU so that the number of possible potential 
energies of the system between U m \ n and C/ max , to a precision of AU, is (t/ max — U m \ n )/ AU 
and becomes polynomial in the "expensive" variables. The smallest possible change in 
potential energy between two different configurations of the particles occurs roughly when 
two electrons are separated by a large distance and one of them is very slightly displaced. 
The maximum distance between the two electrons is of order L, and the minimum dis- 
placement is of order 5. Let the initial r vector between the electrons have length L and 
let the final one have length V . The minimum change in distance V — L occurs when 
the displacement vector of length 5 makes approximately a right angle with the initial r 
vector, so that V = \J L 2 + 5 2 « L + 5 2 /2L. This assumption is justified because the angle 
6 between the initial r vector and the displacement vector for which V = L is such that 
cos# = 5/2L 1, so the angle that minimizes V — L can be made arbitrarily close to ir/2. 
Then the minimum change in potential energy due to this small displacement is 
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whence 

U^-U^ ^ ^!p (iVe) Zl) . . . ? Zjv J = ^3p (jVe) _ _ _ ? (13) 

where P(N e , Zi, . . . , Zjv„) is some polynomial expression of the arguments with constant 
factors absorbed into it. Hence it may be possible to construct the potential energy op- 
erators in a number of gates polynomial in iV e , the Zi, and N as long as we restrict the 
precision of the potential energy values to AU and utilize redundancy effectively. 

3 Numerical Tests 

We illustrate the use of these operators in an actual algorithm with two examples. 
3.1 Order of Error: Particle in a Box 

We first numerically demonstrate the accuracy of our approximation for e~n Te by using 
it to determine a time evolution for a particle in a one-dimensional infinite potential well 
of length L, hence dropping the Coulombic term from the Hamiltonian. The goal here is 
to determine the order of error of our approximations by simulating the simplest possible 
system. We enforce the boundary conditions by introducing a potential energy operator 
that takes an arbitrarily large value at the endpoints. 

An exact wavefunction beginning in a uniform superposition can be expressed as a linear 
combination of eigenfunctions by Fourier expansion: 

I 2 3 / 2 00 1 

^exact(x,0) = —= (0 < X < L) => Vexact(M) = Y\ ^ r V>2fc-1 ^e' 1 ^^ 

k=l 

where ip a = y/2/L sm(airx / L) , with corresponding energies E a = a 2 h 2 7r 2 /2mL 2 . We com- 
pare Inexact t)\ 2 to the probability density of an initial input \ip(x,0)) = ^772 X^^=o 1 K) 
to the quantum algorithm as it evolves under the time evolution operators described in 
Section [2] Figure [2] shows qualitatively the probability density of ^ ex act(^,t) versus that 
of \ip(x,t)) for several evolution times T. In all simulations, 1000 terms are used in the 
Fourier series to compute tp e ^a,ct(x, t). 

To estimate the convergence properties of the stated algorithm, we compute the root 
mean square error of the probability density of \ip(x, t)) produced by the quantum algorithm 
versus both spatial and temporal discretization: 

^1/2 

RMSE = -f-l V \(^x,t)\ij(x,t))-\i; exact (x,t)f^ 




The plots of Figure [3] indicate that the error scales as a power of 5, while preliminary results 
do not reveal the order of error in the temporal discretization. We additionally compute 
the measure of error used by Yepez and Boghosian, defined by E YB = 2" n / 2 (RMSE) ([13]. 
eq. 30). The slope of log(RMSE) in the left plot is 0.2547 while that of log(£ Y B) is 0.7547, 
so the respective errors scale as 0(<5 ' 2547 ) and O(5 7547 ). This convergence in space is 
not nearly as efficient as with the QLGA representation of Schrodinger evolution, which 



8 



T = 1 x 10" 3 T = 2 x 1(T 3 T = 3 x 1(T 3 




0.1 2 0.3 4 OS OB 0.7 8 0.9 1 0.1 2 0.3 4 5 OB 0.7 8 9 1 0.1 2 3 0.4 5 6 7 8 9 1 
2 | 1 1 1 1 1 1 1 , 1 1 2.5 , 1 1 1 1 1 1 1 , 1 1 2.5 , , 1 1 1 , 1 , 1 1 1 




0.1 2 0.3 4 5 6 0.7 08 9 1 0.1 2 0.3 4 5 6 0.7 8 9 1 0.1 2 3 0.4 5 6 7 8 9 



Figure 2: Comparison of probability densities for a particle in a box as predicted by the 
quantum algorithm (top row) and from analytical solutions (bottom row), evolved from a 
"flat" distribution. We use N t = 1000, n = 10, and L = 1. Time is in atomic units. 



usually converges with fourth-order error (Eyb) in <5 [13J. Thus future work might seek to 
increase the efficiency of our intuitive "Cartesian" operator representation. From the right 
plot, the envelope of the maximum RMSE also appears to decrease slightly sublinearly 
with decreasing e, while interesting nonlinear patterns appear that suggest the presence of 
vertical asymptotes and more complicated behavior for most e. Of course, the temporal 
error in the time evolution operator itself is second-order by equation ([5]), although this 
error can be reduced at the cost of increased gate complexity with a third-order Trotter 
decomposition. 

Several subtle sources of numerical error arise in our simulations. One downside to 
our construction of the kinetic energy operator is that its endpoint error, as mentioned in 
Section 2^3 propagates with each application of the time evolution operator constructed 
from the approximate P in equation Q. However, this error affects at most a fixed small 
number of basis states in the particle state vectors at each time step. As long as the spatial 
discretization sufficiently exceeds the temporal discretization, or Nt <C 2 n , it should not 
noticeably hinder the simulation. Figure [2] shows that even when Nt ~ 2 n , one obtains 
qualitatively excellent results: hence our construction of P as diagonal in the momentum 
basis is indeed robust. Conversely, if \h) differs only slightly from an eigenvector of H (as 
a result of the approximation for P) , then e lcH \ h) poorly approximates a phase shift of | h) 
for large c, so short absolute simulation times T confer the best numerical stability. 

Reasons that the RMSE does not converge exactly to when either of the discretizations 
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Figure 3: Left: log-log plot of RMSE of \ip(x, t)\ 2 versus spatial interval 5 at fixed Nt = 1000 
for n = 1,2, . . . , 10. Both log(RMSE) and log(i?YB) are shown as functions of log 5, with 
slopes of 0.2547 and 0.7547, respectively. Right: RMSE versus time step e at fixed n = 6 
for 10 < Nt < 1000. In both plots, T = 10~ 3 and L = 1. The error does not converge to 
at spatial or temporal discretization due to the fixed error in the other discretization. 



tend to include the fixed error in the other variable, endpoint error, error due to a finite 
boundary potential, and (to a much lesser degree) error from a finite Fourier sum. 

3.2 Simple Molecules 

Though the electronic structure of molecules represents a physically interesting case, the 
present algorithm is limited in this domain because modeling electrons and nuclei as spinless 
precludes the simulation of atomic orbitals beyond the Is. 

Nonetheless, the algorithm is simulated in MATLAB to determine what it would predict 
for the electronic charge density of several molecules, each with at most two spinless "elec- 
trons" and three nuclei. To emulate the physical setup of a molecular system, the electron 
wavefunctions are initialized as uniform superpositions over subsets of position basis states 
via Hadamard transforms, while the much more massive nuclei are initialized to single po- 
sition basis states. For the range of precision that we consider, antisymmetrization makes 
little noticeable difference. For simplicity, the numerical implementation uses fully quantum 
encoding and operators but the classical clamped-nucleus Born-Oppenheimer approxima- 
tion. We thus include only the terms T e , U ee , and U en in the molecular Hamiltonian ([!]). 
Also, we carry out the simulation in two dimensions both to reduce the space complexity 
and for easier visualization. A normalized time interval of T = 1 is used. The results of the 
algorithm for four very simple molecules are displayed in Figure |4j These plots show effec- 
tive cross-sections of three-dimensional orbitals, where each square represents an electron 
position basis state. 

The very limited precision of the classical computer does not provide a detailed picture, 
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Figure 4: 2-D electron densities of H, H2 , H2, and H3 , respectively, generated by the 
quantum simulation algorithm. 



but the plots all exhibit symmetry with respect to the nuclei. To simulate H2 and H3 
efficiently, the precision was decreased by a factor of four. In these latter two plots, the 
greatest electron density appears to be in the center, where the atomic orbitals overlap. 
The simulated hydrogen and orbitals, though lobed and artificially resembling those of 
a 3d subshell, most likely reflect numerical errors in the simulation due to discretization. 
Compare the plots of Figure [4] to the classical results of Figure [5j 




Figure 5: Analytical and classical results for \ip\ 2 of H, H^, H2, and (only the latter 
two molecules are in the ground state). The first image shows the exact electron density of 
an H atom corresponding to a 3d orbital in two dimensions, with wavefunction of the form 
ip(ro,0) = r]rQe~ r °/ 3 sin 20. is a linear combination of these exact 3d H orbitals. The 
last two plots are outputs of classical variational methods. 

We have not addressed the question of simulating both ground and excited states of 
systems such as molecules. Adjusting the input qubits to represent an initial configuration 
with a desired energy might allow simulation of states with arbitrary energies. Alternatively, 
Zalka's proposed "energy drain" method could induce a simulated system to decay into its 
ground state by coupling it to an external reservoir [3]. 

4 Conclusion 

We have constructed the relevant quantum gates for simulating the non-relativistic Coulomb 
Hamiltonian for spinless many-body quantum systems, building on Wiesner and Zalka's 
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model. We have suggested implementations of these gates and proven the exactness of the 
kinetic energy operator in the continuum limit. Lastly, we have simulated the algorithm in 
MATLAB and compared its outputs to analytical solutions, finding good agreement. More 
numerical tests are needed and forthcoming. In particular, it would be desirable to simulate 
two particles acting via the Coulomb force in a potential well in 1, 2, and 3 dimensions, 
and to compare those numerical tests involving interacting particles to classically obtained 
benchmarks. 

Our preliminary investigation of this algorithm, under highly simplifying assumptions, 
has two main limitations. First, disregarding fermionic spin statistics (the spin component 
of the wavefunction) means there is no way to enforce the Pauli exclusion principle, which 
does not currently permit the recovery of realistic atomic and molecular structure. Second, 
if one chooses to sample the probability density of the many-body wavefunction at the end 
of the time evolution by making repeated runs of the algorithm and measuring (Section 



2.2), the number of runs necessary to obtain sufficiently good "resolution" is unknown. 
This latter problem is also encountered in the algorithm of Yepez and Boghosian [13] . 

It should be possible to incorporate smaller terms into the Hamiltonian with no change 
in the form of the momentum operator. Doing so lifts many of the simplifications we have 
made. For example, the Dirac Hamiltonian corrects for both spin and relativistic effects 
|18j . Two extra qubits are required to represent positive- and negative-energy states as well 
as spin up or down. An electron wavefunction would thus reside in a Hilbert space Hr^Hf 2 
where 7i r is 2 3n -dimensional and 7i s is two-dimensional. This results in four-component 
spinor wavefunctions where each component is itself a position-dependent wavefunction. 
The Dirac Hamiltonian for an electron in an electromagnetic field is 

H D = 7 °(mc 2 + yV M c) + U 

where 7 M are the gamma matrices and tt = P + eA / c where P and A are operator- valued 
vectors of momentum operators and vector potential operators, respectively, in the x, y, 
and z directions. The momentum operators in each of the three spatial dimensions act only 
on the qubits encoding the electron's position in that dimension. Namely, P = (P x , P y , P z ) 
with P x = P (g> J 2 2n, P y = h n ® P ® h n , Pz = l2 2n ® P-i an d P as defined in equation Q. 
Because the gamma matrices act on the spinor components, we can rewrite these operators 
in terms of matrices of the correct dimension as 

^relativistic + spin = l°mc 2 ® I 2 3n + (7°7 M ® 7 2 3 ")(^4 ® TT/^c) +h®U. 

Both the position and the velocity of every particle in the system make a nonlocal contri- 
bution to the vector potential at every point, so implementing A might therefore require 
additional qubits for keeping track of such velocities. Finally, these corrections we have 
described account only for electron spin and not nuclear spin. 
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